library("boot")

# bootstrap function features
boot.mean <- function(data, i) {
  x <- data[i]
  return(mean(x))
}

# calculate swept area density from input
sweptarea.density <- function(C, S, k, R = 10000) {
  # get density as d = c / s / k
  d = C / S / k
  # bootstrap d
  boot.d <- boot(d, statistic=boot.mean, R = R)
  return(boot.d)
}

# calc CI 
sweptarea.ci <- function(density, level = 0.95) {
  return(boot.ci(density, conf = level, type="all"))
}

#' Swept area method function
#' @param C Catch weight volume
#' @param S Fishing space area square
#' @param k Fishing catchibility coefficient
#' @param total_space Total interpolation water area square
#' @param dim Denumerator for catch weights, eg. if input data in grams 1000000 will lead result to tones
#' @param p Confidence level
#' @param plot Display result fig, TRUE/FALSE
sweptarea <- function(C, S, k, total_space, dim = 1, R = 10000, p = 0.95, plot = FALSE, func = "median") {
  # fit density bootstrap
  d <- sweptarea.density(C = C, S = S, k = k, R = R)
  # get density CI
  ci <- sweptarea.ci(density = d, level = p)
  print(ci)
  est <- NULL
  if (func == 'mean') {
    est <- mean(d$t)
  } else {
    est <- median(d$t)
  }
  
  # estimate mediam biomass from density and total_space area
  b <- est * total_space / dim
  # estimate biomass CI
  b.lo <- ci$bca[4] * total_space / dim
  b.hi <- ci$bca[5] * total_space / dim
  res <- list(B = b, B_CI_LOW = b.lo, B_CI_HIGH = b.hi, d.fit = d, d.ci = ci)
  if (plot) {
    plot(density(d$t * total_space / dim), main = "Biomass normal distribution approximation", xlab = "B, biomass")
    abline(v = b, col="blue", lty="dashed")
    #abline(v = b.lo, col="red", lty = "dashed")
    #abline(v = b.hi, col="red", lty="dashed")
  }
  
  
  return(res)
}
